
model{

#== pop 1 ==#

for(i in 1:N){
for (k in 1:K){
s1_1[i,k]<- se[k]^x1[i,k]*((1-se[k])^(1-x1[i,k]))
s2_1[i,k]<- sp[k]^(1-x1[i,k])*((1-sp[k])^x1[i,k])
}
for (j in 1:K){
for (h in 1:K){
cop1[i,j,h]<- c1[j,h]*(-1)^(x1[i,j] + x1[i,h])/(s1_1[i,j]*s1_1[i,h])
con1[i,j,h]<- c2[j,h]*(-1)^(x1[i,j] + x1[i,h])/(s2_1[i,j]*s2_1[i,h])
}
}
eta1[i] = (prod(s1_1[i,1:K]) *(1+ sum(cop1[i,,])))
theta1[i] =(prod(s2_1[i, 1:K]) *(1+sum(con1[i,,])))
prob1[i]=pi1*eta1[i] + (1-pi1)*theta1[i]
z1[i] ~ dpois( - log(prob1[i]))
}
t1[1:N] ~ dmulti(prob1[1:N], n1)

#== pop 2 ==#

for(i in 1:N){
for (k in 1:K){
s1_2[i,k]<- se[k]^x2[i,k]*((1-se[k])^(1-x2[i,k]))
s2_2[i,k]<- sp[k]^(1-x2[i,k])*((1-sp[k])^x2[i,k])
}
for (j in 1:K){
for (h in 1:K){
cop2[i,j,h]<- c1[j,h]*(-1)^(x2[i,j] + x2[i,h])/(s1_2[i,j]*s1_2[i,h])
con2[i,j,h]<- c2[j,h]*(-1)^(x2[i,j] + x2[i,h])/(s2_2[i,j]*s2_2[i,h])
}
}
eta2[i] = (prod(s1_2[i,1:K]) *(1+ sum(cop2[i,,])))
theta2[i] =(prod(s2_2[i, 1:K]) *(1+sum(con2[i,,])))
prob2[i]=pi2*eta2[i] + (1-pi2)*theta2[i]
z2[i] ~ dpois( - log(prob2[i]))
}
t2[1:N] ~ dmulti(prob2[1:N], n2)


for (k in 1:K){
se[k] ~ dbeta( omega1*(kappa1 -2)+1, (1-omega1)*(kappa1-2) +1)
sp[k] ~ dbeta( omega2*(kappa2 -2)+1, (1-omega2)*(kappa2-2) +1)
}
for (l in 1:(K-1)){
for (h in (l+1):K){
c1[l,h] ~ dunif((se[l]-1)*(1-se[h]), (min(se[l],se[h])-se[l]*se[h]))
c2[l,h] ~ dunif((sp[l]-1)*(1-sp[h]), (min(sp[l],sp[h])-sp[l]*sp[h]))
}}
for (h in 1:K){
for (l in h:K){
c1[l,h] <-0
c2[l,h] <-0
}}
omega1 ~ dbeta(1,1)T(0.5,)
omega2 ~ dbeta(1,1)T(0.5,)
kappa1 = kappaMinusTwo1 +2
kappaMinusTwo1~ dgamma(0.01,0.01)
kappa2 = kappaMinusTwo2+2
kappaMinusTwo2 ~ dgamma(0.01,0.01)
pi1 ~ dbeta(1.8,2)
pi2 ~ dbeta(1.8,2)
}
